
#include <stdio.h>
#include <inttypes.h>

double invSqrtd(double x) 
{
  double halfx = 0.5 * x;
  double y = x;
  uint64_t i = *(uint64_t*)&y;
  i = 0x5fe6ec85e7de30da - (i>>1); // the magic number
  y = *(double*)&i;
  y *= (1.5 - (halfx * y * y));    // 1st iteration
  y *= (1.5 - (halfx * y * y));    // 2nd iteration
  return y;
}

float invSqrtf(float x) 
{
	float halfx = 0.5f * x;
	float y = x;
	long i = *(long*)&y;
	i = 0x5f3759df - (i>>1);
	y = *(float*)&i;
	y = y * (1.5f - (halfx * y * y));
	return y;
}

int main(int argc, char* argv[])
{
  double a = 16.0;
  float af = 16.0f;
  float bf = invSqrtf(af);
  double b = invSqrtd(a);

  printf("%f %lf\n", bf, b);

  return 0;
}
